/* Appliation of velocity boundary condition. */
# include "scene.h"
# include <stdio.h>
# include <stdlib.h>
# include <math.h>

# ifdef _OMP_VEL_BC_
# include <omp.h>
# endif

extern double inlet_velocity(double r, double t, short index);

extern void vel_bc(struct grid *g, short tflag, unsigned short iRK)
  {
   short nr, nth, nz, i, _i, j, _j, k;
   double ***u, ***v, ***w, *r, *z, *th, t;

#  ifdef TRACE_FILE
   fprintf(g->trptr, "Entering into 'vel_bc()'...\n");
#  endif

   if( ((iRK > 0)&&(tflag != 0)) || ((iRK == 0)&&(tflag != 1)) || (iRK > 2) )
       {
        fprintf(g->lptr, "\nERROR: In 'vel_bc()': (tflag, iRK) = (%d, %d).\nThis function is designed to handle only following cases:\n\t[tflag = 0, iRK = 1, 2] and [tflag = 1, iRK = 0]\n", tflag, iRK);
        exit(-1);
       }

   t = g->t + tflag*g->dt;

   u = g->u;   v = g->v;   w = g->w;   nr = g->nr;   nth = g->nth;   nz = g->nz;
   r = g->r;   th = g->th; z = g->z;

/*** w[i][j][NzG-2] is needed for application of the Numan Pressure-BC at the bottom interface and also to
   evaluare RHS of the PPE. ***/
# ifdef _OMP_VEL_BC_
  #pragma omp parallel private(i,j,k,_i,_j)
# endif
    {
  /* At the TOP-SURFACE */
#    ifdef _OMP_VEL_BC_
     #pragma omp for schedule(guided) nowait
#    endif
      for(i = NrG; i < nr+NrG; i++) for(j = NthG; j < NthG+nth; j++)
         {
          u[i][j][NzG+nz] = -u[i][j][NzG+nz-1];
          v[i][j][NzG+nz] = -v[i][j][NzG+nz-1];
          w[i][j][NzG+nz] = w[i][j][NzG+nz-2];
         }

 /*  At the SIDE-WALL */
#    ifdef _OMP_VEL_BC_
     #pragma omp for schedule(guided) nowait
#    endif
     for(k = NzG; k < NzG + g->nz; k++) for(j = NthG; j < NthG+nth; j++)
         {
          u[NrG+nr][j][k] = (r[NrG+nr-1]/r[NrG+nr+1])*u[NrG+nr-2][j][k];
          v[NrG+nr][j][k] = -v[NrG+nr-1][j][k];
          w[NrG+nr][j][k] = -w[NrG+nr-1][j][k];
         }

  /* At the INLET */
#    ifdef _OMP_VEL_BC_
     #pragma omp for schedule(guided) nowait
#    endif
     for(i = NrG; i < NrG + g->nr_nozzle; i++) for(j = NthG; j < NthG+nth; j++)
         {
          u[i][j][NzG-1] = -u[i][j][NzG];
          v[i][j][NzG-1] = -v[i][j][NzG];

          switch(iRK)
               {
                case 0: w[i][j][NzG-1] = inlet_velocity(0.5*(r[i]+r[i+1]), t, 0);
                        break;

                case 1: w[i][j][NzG-1] = inlet_velocity(0.5*(r[i]+r[i+1]), t, 0)\
                             + c2*g->dt*inlet_velocity(0.5*(r[i]+r[i+1]), t, 1);
                        break;

                case 2: w[i][j][NzG-1] = inlet_velocity(0.5*(r[i]+r[i+1]), t, 0)\
                             + (a31+a32)*g->dt*inlet_velocity(0.5*(r[i]+r[i+1]), t, 1)\
                             + a32*c2*g->dt*g->dt*inlet_velocity(0.5*(r[i]+r[i+1]), t, 2);
                        break;
               }

          w[i][j][NzG-2] = w[i][j][NzG];
         }

  /* At the SOLID-BASE after the inlet and before the outlet. */
#    ifdef _OMP_VEL_BC_
     #pragma omp for schedule(guided) nowait
#    endif
      for(i = NrG+g->nr_nozzle; i < NrG+nr-g->nr_gap; i++) for(j = NthG; j < NthG+nth; j++)
         {
          u[i][j][NzG-1] = -u[i][j][NzG];
          v[i][j][NzG-1] = -v[i][j][NzG];
          w[i][j][NzG-1] = 0;
          w[i][j][NzG-2] = w[i][j][NzG];
         }

  /* At the OUTLET: Orlansky BC */
#    ifdef _OMP_VEL_BC_
     #pragma omp for schedule(guided)
#    endif
     for(_i = 0; _i < g->nr_gap; _i++) for(i = _i+NrG+nr-g->nr_gap, j = NthG, _j = 0; _j < nth; _j++, j++)
         {
          switch(iRK)
            {
             case 0:
             u[i][j][NzG-1] = g->outBC.tnext_ug[_i][_j];
             v[i][j][NzG-1] = g->outBC.tnext_vg[_i][_j];
             w[i][j][NzG-1] = g->outBC.tnext_wg[_i][_j];
             break;

             case 1:
             u[i][j][NzG-1] = c2*(g->outBC.u[2][_i][_j][0]/6 - g->outBC.u[1][_i][_j][0] + g->outBC.u[0][_i][_j][0]/2\
                              + g->outBC.tnext_ug[_i][_j]/3) + g->outBC.u[0][_i][_j][0];

             v[i][j][NzG-1] = c2*(g->outBC.v[2][_i][_j][0]/6 - g->outBC.v[1][_i][_j][0] + g->outBC.v[0][_i][_j][0]/2\
                              + g->outBC.tnext_vg[_i][_j]/3) + g->outBC.v[0][_i][_j][0];

             w[i][j][NzG-1] = c2*(g->outBC.w[2][_i][_j][0]/6 - g->outBC.w[1][_i][_j][0] + g->outBC.w[0][_i][_j][0]/2\
                              + g->outBC.tnext_wg[_i][_j]/3) + g->outBC.w[0][_i][_j][0];
             break;

             case 2:
             u[i][j][NzG-1] = (a31+a32)*g->outBC.u[2][_i][_j][0]/6 - (a31+a32-a32*c2)*g->outBC.u[1][_i][_j][0]\
                 + (1 + 0.5*(a31+a32-4*a32*c2))*g->outBC.u[0][_i][_j][0] + (a31+a32+3*a32*c2)*g->outBC.tnext_ug[_i][_j]/3;

             v[i][j][NzG-1] = (a31+a32)*g->outBC.v[2][_i][_j][0]/6 - (a31+a32-a32*c2)*g->outBC.v[1][_i][_j][0]\
                 + (1 + 0.5*(a31+a32-4*a32*c2))*g->outBC.v[0][_i][_j][0] + (a31+a32+3*a32*c2)*g->outBC.tnext_vg[_i][_j]/3;

             w[i][j][NzG-1] = (a31+a32)*g->outBC.w[2][_i][_j][0]/6 - (a31+a32-a32*c2)*g->outBC.w[1][_i][_j][0]\
                 + (1 + 0.5*(a31+a32-4*a32*c2))*g->outBC.w[0][_i][_j][0] + (a31+a32+3*a32*c2)*g->outBC.tnext_wg[_i][_j]/3;

             break;
            }
         }

  /* At the corner-cell after the outlet. */
#    ifdef _OMP_VEL_BC_
     #pragma omp for schedule(guided) nowait
#    endif
/*** This requires the outlet-BC to be applied first. ***/
     for(j = NthG; j < NthG+nth; j++) w[NrG+nr][j][NzG-1] = -w[NrG+nr-1][j][NzG-1];

  /* For the NEGATIVE radius: applying symmetric B.C. across the center. */
#    ifdef _OMP_VEL_BC_
     #pragma omp for schedule(guided)
#    endif
/*** The i-variation should always be kept in the end. ***/
     for(k = 0; k < 2*NzG+nz; k++) for(j = NthG; j < NthG+nth; j++) for(i = NrG; i >= 1; i--)
         {
     /*** Order of this if-else should not be reversed. ***/
          u[NrG-i][j][k] = (i > 1 ? -u[NrG+i-2][j][k]:0.5*(u[NrG-2][j][k]+u[NrG][j][k]));
          v[NrG-i][j][k] = -v[NrG+i-1][j][k];
          w[NrG-i][j][k] = w[NrG+i-1][j][k];
         }


/*   For ghost cells in theta-direction. */
#    ifdef _OMP_VEL_BC_
     #pragma omp for schedule(guided)
#    endif
  /* It needs the other boundary conditions applied first. */
     for(i = 0; i < nr+2*NrG; i++) for(k = 0; k < nz+2*NzG; k++) for(j = 0; j < NthG; j++)
        {
      /* Applying symmetric B.C. */
         u[i][NthG-1-j][k] = u[i][NthG+j][k];
         u[i][NthG+nth+j][k] = u[i][NthG+nth-1-j][k];

         v[i][NthG-1-j][k] = v[i][NthG+j][k];
         v[i][j+NthG+nth][k] = v[i][NthG+nth-1-j][k];

         w[i][NthG-1-j][k] = w[i][NthG+j][k];
         w[i][j+NthG+nth][k] = w[i][NthG+nth-1-j][k];
        }

/*** Although we don't apply Numan Pressure-BC at the OUTLET, we need w[i][j][NzG-2] because it is needed to
     evaluate the R.H.S. of PPE. It is evaluated from the divergence-free condition of the cells at the outlet
     boundary. It needs other BCs applied first, inluding the BC in theta-direction. ***/
#    ifdef _OMP_VEL_BC_
     #pragma omp for schedule(guided)
#    endif
     for(i = NrG+nr-g->nr_gap; i < NrG+nr; i++) for(j = NthG; j < NthG+nth; j++)
     w[i][j][NzG-2] = w[i][j][NzG] + (z[NzG+1]-z[NzG-1])*(\
                        (r[i+1]*(u[i][j][NzG]+u[i][j][NzG-1]) - r[i]*(u[i-1][j][NzG]+u[i-1][j][NzG-1]))/(r[i+1]-r[i])\
                      + (v[i][j][NzG]+v[i][j][NzG-1]-v[i][j-1][NzG]-v[i][j-1][NzG-1])/(th[j+1]-th[j])\
                                                           )/(r[i]+r[i+1]);

    }

#  ifdef TRACE_FILE
   fprintf(g->trptr, "...'vel_bc()' ends.\n");
#  endif

   return;
  }
